Offline Data Analysis of Strongyloides venezuelensis bulk RNAseq data

Introduction

The goal of this file is to lay out proof-of-concept analyses of the Strongyloides papillosus RNAseq dataset originally analyzed by Hunt et al 2016.

This dataset includes samples representing 5 life stages, each with 2 biological replicates: L1/L2s, FLF, FLM, iL3, PF.

Data Pre-Processing

A full description of Kallisto alignment and data filtering/normalization steps can be found in Sp_RNAseq_Data_Processing.rmd.

Data Analysis

The limma package ( Ritchie et al 2015, Phipson et al 2016) is used to conduct pairwise differential gene expression analyses between life stages. The results of the pairwise comparison is displayed as a volcano plot and interactive DataTable.

Code

All code is echoed under descriptive headers; code chunks are hidden from view by default. Users may show hidden R code by clicking the Show buttons. In addition, all code chunks are collated at the end of the document in an Appendix.

Load and Parse Preprocessed Data

This code loads R data objects that has been preproprocessed by Sp_RNAseq_Data_Processing.rmd.

# Load and Parse Preprocessed Data
load (file = "../Outputs/SpRNAseq_data_preprocessed")
targets <- SpRNAseq.preprocessed.data$targets
annotations <- SpRNAseq.preprocessed.data$annotations
log2.cpm.filtered.norm <- SpRNAseq.preprocessed.data$log2.cpm.filtered.norm
myDGEList.filtered.norm <-SpRNAseq.preprocessed.data$myDGEList.filtered.norm

Loading required package: edgeR
Loading required package: limma

rm(SpRNAseq.preprocessed.data)

load(file = "../Strongyloides_RNAseq_Browser/Data/Sp_vDEGList")

Hierarchical Clustering and Principle Components Analysis

This code chunk starts with filtered and normalized abundance data in a data frame (not tidy). It will implement hierarchical clustering and PCA analyses on the data. It will plot various graphs, including a dendrogram of the heirachical clustering, and several plots of visualize the PCA. Because the data that are passed into these analyses do not have batch correction applied, the clustering appears dominanted by a batch effect.

# Introduction to this chunk -----------
# This code chunk starts with filtered and normalized abundance data in a data frame (not tidy).
# It will implement hierarchical clustering and PCA analyses on the data.
# It will plot various graphs and can save them in PDF files.
# Load packages ------
suppressPackageStartupMessages({
  library(tidyverse) # you're familiar with this fromt the past two lectures
  library(ggplot2)
  library(RColorBrewer)
  library(ggdendro)
  library(magrittr)
  library(factoextra)
  library(gridExtra)
  library(cowplot)
  library(dendextend)
})

# Identify variables of interest in study design file ----
group <- factor(targets$group)

# Hierarchical clustering ---------------
# Remember: hierarchical clustering can only work on a data matrix, not a data frame

# Calculate distance matrix
# dist calculates distance between rows, so transpose data so that we get distance between samples.
# how similar are samples from each other
colnames(log2.cpm.filtered.norm)<-paste(targets$group,substr(targets$sample, 9,10), sep = ".")
distance <- dist(t(log2.cpm.filtered.norm), method = "maximum") #other distance methods are "euclidean", maximum", "manhattan", "canberra", "binary" or "minkowski"

# Calculate clusters to visualize differences. This is the hierarchical clustering.
# The methods here include: single (i.e. "friends-of-friends"), complete (i.e. complete linkage), and average (i.e. UPGMA). Here's a comparison of different types: https://en.wikipedia.org/wiki/UPGMA#Comparison_with_other_linkages
clusters <- hclust(distance, method = "complete") #other agglomeration methods are "ward.D", "ward.D2", "single", "complete", "average", "mcquitty", "median", or "centroid"
dend <- as.dendrogram(clusters) 

p1<-dend %>% 
  dendextend::set("branches_k_color", k = 6) %>% 
  dendextend::set("hang_leaves", c(0.05)) %>% 
  dendextend::set("labels_cex", c(0.5)) %>%
  dendextend::set("labels_colors", k = 6) %>% 
  dendextend::set("branches_lwd", c(0.7)) %>% 
  
  as.ggdend %>%
  ggplot (offset_labels = -0.2) +
  theme_dendro() +
  ylim(0, max(get_branches_heights(dend))) +
  labs(title = "S. papillosus: Hierarchical Cluster Dendrogram",
       subtitle = "filtered, TMM normalized",
       y = "Distance",
       x = "Life stage") +
  coord_fixed(1/2) +
  theme(axis.title.x = element_text(color = "black"),
        axis.title.y = element_text(angle = 90),
        axis.text.y = element_text(angle = 0),
        axis.line.y = element_line(color = "black"),
        axis.ticks.y = element_line(color = "black"),
        axis.ticks.length.y = unit(2, "mm"))

Dendrogram to visualize heirachical clustering

Clustering performed on filtered and normalized abundance data using the “complete” method.

# Principal component analysis (PCA) -------------
# this also works on a data matrix, not a data frame
pca.res <- prcomp(t(log2.cpm.filtered.norm), scale.=F, retx=T)
#summary(pca.res) # Prints variance summary for all principal components.

#pca.res$rotation #$rotation shows you how much each gene influenced each PC (called 'scores')
#pca.res$x # 'x' shows you how much each sample influenced each PC (called 'loadings')
#note that these have a magnitude and a direction (this is the basis for making a PCA plot)
## This generates a screeplot: a standard way to view eigenvalues for each PCA. Shows the proportion of variance accounted for by each PC. Plotting only the first 10 dimensions.
p2<-fviz_eig(pca.res,
             barcolor = brewer.pal(8,"Pastel2")[8],
             barfill = brewer.pal(8,"Pastel2")[8],
             linecolor = "black",
             main = "S. papillosus Scree plot: proportion of variance accounted for by each principal component",
             ggtheme = theme_bw()) 

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Screeplot of PCA Eigenvalues

A scree plot is a standard way to view eigenvalues for each PCA. The plot shows the proportion of variance accounted for by each PC.

pc.var<-pca.res$sdev^2 # sdev^2 captures these eigenvalues from the PCA result
pc.per<-round(pc.var/sum(pc.var)*100, 1) # we can then use these eigenvalues to calculate the percentage variance explained by each PC

# Visualize the PCA result ------------------
#lets first plot any two PCs against each other
#We know how much each sample contributes to each PC (loadings), so let's plot
pca.res.df <- as_tibble(pca.res$x)

# Plotting PC1 and PC2
p3<-ggplot(pca.res.df) +
  aes(x=PC1, y=PC2, label=targets$group, 
      fill = targets$group,
      color = targets$group
  ) +
  geom_point(size=4, shape= 21, color = "black", alpha = 0.5) +
  #geom_label(color = "black", size = 2) +
  #scale_fill_brewer(palette = "Set2") +
  #scale_color_brewer(palette = "Set2", guide = FALSE) +
  #stat_ellipse() +
  xlab(paste0("PC1 (",pc.per[1],"%",")")) + 
  ylab(paste0("PC2 (",pc.per[2],"%",")")) +
  labs(title="S. papillosus: Principal Components Analysis of RNAseq Samples",
       caption = "Note: analysis is blind to life stage identity.",
       subtitle ="All samples/groups",
       fill = "Life Stage") +
  scale_x_continuous(expand = c(.3, .3)) +
  scale_y_continuous(expand = c(.3, .3)) +
  coord_fixed() +
  theme_bw()

PCA Plot

Plot of the samples in PCA space. Fill color indicates life stage.

# Create a PCA 'small multiples' chart ----
pca.res.df <- pca.res$x[,1:3] %>% 
  as_tibble() %>%
  add_column(sample = targets$sample,
             group = targets$group)

pca.pivot <- pivot_longer(pca.res.df, # dataframe to be pivoted
                          cols = PC1:PC3, # column names to be stored as a SINGLE variable
                          names_to = "PC", # name of that new variable (column)
                          values_to = "loadings") # name of new variable (column) storing all the values (data)
PC1<-subset(pca.pivot, PC == "PC1")
PC2 <-subset(pca.pivot, PC == "PC2")
#PC3 <- subset(pca.pivot, PC == "PC3")
#PC4 <- subset(pca.pivot, PC == "PC4")

p6<-ggplot(pca.pivot) +
  aes(x=sample, y=loadings) + # you could iteratively 'paint' different covariates onto this plot using the 'fill' aes
  geom_bar(stat="identity", fill = brewer.pal(8,"Pastel2")[8]) +
  scale_fill_brewer(palette = "Set2") +
  facet_wrap(~PC) +
  geom_bar(data = PC1, stat = "identity", aes(fill = group)) +
  geom_bar(data = PC2, stat = "identity", aes(fill = group)) +
  labs(title="S. papillosus: PCA 'small multiples' plot",
       fill = "Life Stage/Source",
       subtitle ="All samples/groups") +
  scale_x_discrete(limits = targets$sample, labels = targets$group) +
  theme_bw() +
  coord_flip()

PCA “Small Multiples” Plot

Heatmap of Gene Expression Across Life Stages

Make a heatmap for all the genes using the Log2CPM values.

suppressPackageStartupMessages({
library(pheatmap)
library(RColorBrewer)
library(heatmaply)
})
# Make a heatmap for all the genes using the Log2CPM values

diffGenes <- v.DEGList.filtered.norm$E %>%
  as_tibble(rownames = "geneID", .name_repair = "unique") %>%
  dplyr::select(!geneID) %>%
  as.matrix()
rownames(diffGenes) <- rownames(v.DEGList.filtered.norm$E)
colnames(diffGenes) <- as.character(v.DEGList.filtered.norm$targets$group)
clustColumns <- hclust(as.dist(1-cor(diffGenes, method="spearman")), method="complete")
clustRows <- hclust(as.dist(1-cor(t(diffGenes), 
                                  method="pearson")), 
                    method="complete") 
par(cex.main=1.2)

showticklabels <- c(TRUE,FALSE)
p<-pheatmap(diffGenes,
            color = RdBu(75),
            cluster_rows = clustRows,
            cluster_cols = clustColumns,
            show_rownames = F,
            scale = "row",
            angle_col = 45,
            main = "S. papillosus: Log2 Counts Per Million (CPM) Expression Across Life Stages"
          
)

Differentially Expressed Genes

This chunk uses a variance-stabilized DGEList of filtered and normalized abundance data. These data/results are examples, a responsive version of this code is avaliable in a Shiny App.

# Introduction to this chunk ----
# Because we have access to biological and technical replicates, we can use statistical tools for differential expression analysis
# Useful reading on differential expression: https://ucdavis-bioinformatics-training.github.io/2018-June-RNA-Seq-Workshop/thursday/DE.html

# Load packages ----
suppressPackageStartupMessages({
  library(tidyverse)
  library(limma) # differential gene expression using linear modeling
  library(edgeR)
  library(gt) 
  library(DT) 
  library(plotly)
  library(ggthemes)
  library(RColorBrewer)
  source("../Strongyloides_RNAseq_Browser/Server/theme_Publication.R")
})

diffGenes.df <- v.DEGList.filtered.norm$E %>%
  as_tibble(rownames = "geneID", .name_repair = "unique")

# Set Expression threshold values for plotting and saving DEGs ----
adj.P.thresh <- 0.05
lfc.thresh <- 1 

group <- factor(v.DEGList.filtered.norm$targets$group)
design <- model.matrix(~0 + group) # no intercept/blocking for matrix, comparisons across group
colnames(design) <- levels(group)


# Fit a linear model to the data ----
fit <- lmFit(v.DEGList.filtered.norm, design = design)

# As an example, generate comparison matrix for a pairwise comparison ----
# iL3s vs FLF
# Note that the target/contrast goups will be divided by the number of life 
# stage groups e.g. PF+FLF/2 - iL3+iL3a+pfL1+ppL1+ppL3/5
comparison <- c('(iL3)-(FLF)')

targetStage<- comparison %>%
  str_split(pattern="-", simplify = T) %>%
  .[,1] %>%
  gsub("(", "", ., fixed = TRUE) %>%
  gsub(")", "", ., fixed = TRUE) %>%
  str_split(pattern = "\\+", simplify = T)

contrastStage<-comparison %>%
  str_split(pattern="-", simplify = T) %>%
  .[,2] %>%
  gsub("(", "", ., fixed = TRUE) %>%
  gsub(")", "", ., fixed = TRUE)  %>%
  str_split(pattern = "\\+", simplify = T)

comparison<- sapply(seq_along(comparison),function(x){
  tS <- as.vector(targetStage[x,]) %>%
    .[. != ""] 
  cS <- as.vector(contrastStage[x,]) %>%
    .[. != ""] 
  paste(paste0(tS, 
               collapse = "+") %>%
          paste0("(",.,")/",length(tS)),
        paste0(cS, 
               collapse = "+") %>%
          paste0("(",.,")/",length(cS)),
        sep = "-")
  
})

# Generate contrast matrix ----
contrast.matrix <- makeContrasts(contrasts = comparison,
                                 levels=design)

# extract the linear model fit -----
fits <- contrasts.fit(fit, contrast.matrix)
# empirical bayes smoothing of gene-wise standard deviations provides increased power (see: https://www.degruyter.com/doi/10.2202/1544-6115.1027)
ebFit <- eBayes(fits)

# Pull out the DEGs that pass a specific threshold for all pairwise comparisons ----
# Adjust for multiple comparisons using method = global. 
results <- decideTests(ebFit, method="global", adjust.method="BH", p.value = adj.P.thresh)

recode01<- function(x){
  case_when(x == 1 ~ "Up",
            x == -1 ~ "Down",
            x == 0 ~ "NotSig")
}
diffDesc <- results %>%
  as_tibble(rownames = "geneID") %>%
  dplyr::mutate(across(-geneID, unclass)) %>%
  dplyr::mutate(across(where(is.double), recode01))

# Function that identifies top DEGs between a specific contrast ----
calc_DEG_tbl <- function (ebFit, coef) {
  myTopHits.df <- limma::topTable(ebFit, adjust ="BH", 
                                  coef=coef, number=40000, 
                                  sort.by="logFC") %>%
    as_tibble(rownames = "geneID") %>%
    dplyr::rename(tStatistic = t, LogOdds = B, BH.adj.P.Val = adj.P.Val) %>%
    dplyr::relocate(UniProtKB, Description, InterPro, GO_term, 
                    In.subclade_geneID, In.subclade_percent_homology,
                    Out.subclade_geneID, Out.subclade_percent_homology,
                    Ce_geneID, Ce_percent_homology, .after = LogOdds)
  
  myTopHits.df
}

list.myTopHits.df <- sapply(comparison, function(y){
  calc_DEG_tbl(ebFit, y)}, 
  simplify = FALSE, 
  USE.NAMES = TRUE)

list.myTopHits.df <- sapply(comparison, function(y){
  list.myTopHits.df[[y]] %>%
    dplyr::select(geneID, 
                  logFC, 
                  BH.adj.P.Val:Ce_percent_homology)},
  simplify = FALSE, 
  USE.NAMES = TRUE)

# Get log2CPM values and threshold information for genes of interest
list.myTopHits.df <- sapply(seq_along(comparison), function(y){
  tS<- targetStage[y,][targetStage[y,]!=""]
  cS<- contrastStage[y,][contrastStage[y,]!=""]
  
  concat_name <- function(x) {
    ifelse(x == "target", 
           paste(tS, collapse = "+"), 
           paste(cS, collapse = "+"))
  }
  
  groupAvgs <- diffGenes.df %>%
    dplyr::select(geneID, starts_with(paste0(tS,"-")), 
                  starts_with(paste0(cS,"-"))) %>%
    pivot_longer(cols = -geneID, names_to = c("group","sample"), values_to = "CPM",
                 names_sep = "-") %>%
    dplyr::mutate(contrastID = if_else(group %in% tS,"target", "contrast")) %>%
    group_by(geneID, contrastID) %>%
    dplyr::select(-sample) %>%
    summarize(mean = mean(CPM), .groups = "drop_last") %>%
    pivot_wider(names_from = contrastID, values_from = mean) %>%
    dplyr::relocate(contrast, .after = target) %>%
    dplyr::rename_with(concat_name, -geneID) %>%
    dplyr::rename_with(.cols =-geneID, .fn = ~ paste0("avg_(",.x,")"))
  
  diffGenes.df %>%
    dplyr::select(geneID, starts_with(paste0(tS,"-")), 
                  starts_with(paste0(cS,"-"))) %>%
    left_join(groupAvgs, by = "geneID") %>%
    left_join(list.myTopHits.df[[y]],., by = "geneID") %>%
    left_join(dplyr::select(diffDesc,geneID,comparison[y]), by = "geneID") %>%
    dplyr::rename(DEG_Desc=comparison[y]) %>%
    dplyr::relocate(DEG_Desc) %>%
    dplyr::relocate(logFC:Ce_percent_homology, .after = last_col())
  
},
simplify = FALSE)

comparison <- gsub("/[0-9]*","", comparison)
names(list.myTopHits.df) <- comparison

list.myTopHits.df <- sapply(comparison, function(y){
  list.myTopHits.df[[y]] %>%
    dplyr::mutate(DEG_Desc = case_when(DEG_Desc == "Up" ~ paste0("Up in ", str_split(y,'-',simplify = T)[1,1]),
                                       DEG_Desc == "Down" ~ paste0("Down in ", str_split(y,'-',simplify = T)[1,1]),
                                       DEG_Desc == "NotSig" ~ "NotSig")) 
},
simplify = FALSE, 
USE.NAMES = TRUE)

# PC1 Volcano Plot and Interactive Table ----
vplot1 <- ggplot(list.myTopHits.df[[1]]) +
  aes(y=-log10(BH.adj.P.Val), x=logFC, text = paste(geneID, "<br>",
                                                    "logFC:", round(logFC, digits = 2), "<br>",
                                                    "p-val:", format(BH.adj.P.Val, digits = 3, scientific = TRUE))) +
  geom_point(size=2) +
  geom_hline(yintercept = -log10(adj.P.thresh), 
             linetype="longdash", 
             colour="grey", 
             size=1) + 
  geom_vline(xintercept = lfc.thresh, 
             linetype="longdash", 
             colour="#BE684D", 
             size=1) +
  geom_vline(xintercept = -lfc.thresh, 
             linetype="longdash", 
             colour="#2C467A", 
             size=1) +
  labs(title = paste0('S. papillosus, Pairwise Comparison: ',
                      gsub('-',
                           ' vs ',
                           comparison[1])),
       subtitle = paste0("grey line: p = ",
                         adj.P.thresh, "; colored lines: log-fold change = ", lfc.thresh),
       color = "GeneIDs") +
  theme_Publication() 
vplot1

# Interactive Tables
yy<- 1
tS<- targetStage[yy,][targetStage[yy,]!=""]
cS<- contrastStage[yy,][contrastStage[yy,]!=""]
sample.num.tS <- sapply(tS, function(x) {colSums(v.DEGList.filtered.norm$design)[[x]]}) %>% sum()
sample.num.cS <- sapply(cS, function(x) {colSums(v.DEGList.filtered.norm$design)[[x]]}) %>% sum()


n_num_cols <- sample.num.tS + sample.num.cS + 5
index_homologs <- length(colnames(list.myTopHits.df[[yy]])) - 5


LS.datatable <- list.myTopHits.df[[yy]] %>%
  DT::datatable(rownames = FALSE,
                caption = htmltools::tags$caption(
                  style = 'caption-side: top; text-align: left; color: black',
                  htmltools::tags$b('Differentially Expressed Genes in', 
                                    htmltools::tags$em('S. papillosus'), 
                                    gsub('-',' vs ',comparison[yy])),
                  htmltools::tags$br(),
                  "Threshold: p < ",
                  adj.P.thresh, "; log-fold change > ",
                  lfc.thresh,
                  htmltools::tags$br(),
                  'Values = log2 counts per million'),
                options = list(autoWidth = TRUE,
                               scrollX = TRUE,
                               scrollY = '300px',
                               scrollCollapse = TRUE,
                               order = list(n_num_cols-1, 
                                            'desc'),
                               searchHighlight = TRUE, 
                               pageLength = 25, 
                               lengthMenu = c("5",
                                              "10",
                                              "25",
                                              "50",
                                              "100"),
                               columnDefs = list(
                                 list(
                                   targets = ((n_num_cols + 
                                                 1)),
                                   render = JS(
                                     "function(data, row) {",
                                     "data.toExponential(1);",
                                     "}")
                                 ),
                                 list(
                                   targets = ((n_num_cols + 
                                                 4):(n_num_cols + 
                                                       5)),
                                   render = JS(
                                     "function(data, type, row, meta) {",
                                     "return type === 'display' && data.length > 20 ?",
                                     "'<span title=\"' + data + '\">' + data.substr(0, 20) + '...</span>' : data;",
                                     "}")
                                 ),
                                 list(targets = "_all",
                                      class="dt-right")
                               ),
                               rowCallback = JS(c(
                                 "function(row, data){",
                                 "  for(var i=0; i<data.length; i++){",
                                 "    if(data[i] === null){",
                                 "      $('td:eq('+i+')', row).html('NA')",
                                 "        .css({'color': 'rgb(151,151,151)', 'font-style': 'italic'});",
                                 "    }",
                                 "  }",
                                 "}"  
                               ))
                               
                )) 
LS.datatable <- LS.datatable %>%
  DT::formatRound(columns=c(3:n_num_cols), 
                  digits=3)

LS.datatable <- LS.datatable %>%
  DT::formatRound(columns=c(n_num_cols+2, 
                            index_homologs+1,
                            index_homologs+3), 
                  digits=2)

LS.datatable <- LS.datatable %>%
  DT::formatSignif(columns=c(n_num_cols+1), 
                   digits=3)

LS.datatable

NA

Benchmarking

Tables with data suitable for benchmarking are not currently available.

Functional Enrichment Analysis

This code prerform GSEA using the clusterProfiler library. Ability to do this depends on the availability of gene sets. Major databases (e.g. msigdb don’t seem to have Strongyloides information. They do have C. elegans gene sets, but I’m not convinced the homology information is good enough for the comparison to be unbiased/meaningful. In Hunt et al 2016, there is an Ensembl Compara protein family set; we will use this as the basis for our gene set libraries.
Note that this uses specific transcript information, which I throw out (e.g. SSTP_0001137400.2 is recoded as SSTP_0001137400).

Given a priori defined set of gene S (e.g., genes shareing the same DO category), the goal of GSEA is to determine whether the members of S are randomly distributed throughout the ranked gene list (L) or primarily found at the top or bottom.
There are three key elements of the GSEA method:
Calculation of an Enrichment Score.
The enrichment score (ES) represent the degree to which a set S is over-represented at the top or bottom of the ranked list L. The score is calculated by walking down the list L, increasing a running-sum statistic when we encounter a gene in S and decreasing when it is not. The magnitude of the increment depends on the gene statistics (e.g., correlation of the gene with phenotype). The ES is the maximum deviation from zero encountered in the random walk; it corresponds to a weighted Kolmogorov-Smirnov-like statistic (Subramanian et al. 2005).
Esimation of Significance Level of ES.
The p-value of the ES is calculated using permutation test. Specifically, we permute the gene labels of the gene list L and recompute the ES of the gene set for the permutated data, which generate a null distribution for the ES. The p-value of the observed ES is then calculated relative to this null distribution.
Adjustment for Multiple Hypothesis Testing.
When the entire gene sets were evaluated, DOSE adjust the estimated significance level to account for multiple hypothesis testing and also q-values were calculated for FDR control.

Appendix I: All code for this report